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ABSTRACT 


A history of uniform pseudorandom number generation is 
presented from the experiments of W.S. Gossett up to the 
present day. A new concept of digitized uniform pseudorandom 
number generator (due to K.V. Leung) is described. Various 
models developed from the original concept have been imple- 
mented and tested for comparison with existing systems. 

These tests produce results showing that the new generators 
appear to have some statistical merit. One practical appli- 


cation of uniform pseudorandom number generation is given. 
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CHAPTER 1 


HISTORY OF RANDOM NUMBER GENERATION 


1.1 Pre-Computer 


The generation of sequences of random numbers from the 
uniform distribution is a problem with a great historical 
background. Gossett (1908) was probably the first person 
to actually use random numbers — his general method con- 
sisted of taking a correlation table of heights and left 
middle finger measurements of 3000 criminals (see MacDonell 
(1901)). The table entries were written on 3000 pieces of 
cardboard, shuffled, and drawn at random (with replacement) 
four at a time. The resulting sequence of random numbers 
was used in formulating what is now known as Student's t 
distribution. However, his "random number generator" was 
extremely slow, and it was hard to say when the 3000 pieces 
of cardboard had been well shuffled. Tippett (1927) re- 
placed this system by a table of 10000 four-digit numbers, 
which had been formed by taking 40000 digits "at random" 
from census tables, and grouping these in fours to obtain 
the 10000 numbers. Fisher and Yates (1938) improved on this 


by producing a table of 15000 random sampling numbers. 


Kendall and Babington-Smith (1938) were the first to 
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successfully generate random digits by mechanical means. 
They used a disc which was rotated by an electric motor at 
aecapids rate an a darkened room. This direc was divided into 
ten equal sections on which appeared the digits 0-9 inclusive. 
An electric spark illuminated the disc instantaneously, so as 
to make the disc appear stationary for a moment. When the 
spark occurred, a number was selected from the disc by 

means of a fixed pointer. The intervals between sparks were 
varied by means of a neon lamp in parallel with a condensor 
in an independent electrical circuit. A key tapped by an 
Observer broke the circuil. intermittantly at irregular 
aWCerVvals;, tO EUrther add to the randomness. “A table of 
random numbers generated by this mechanism was published in 
1939 — this satisfied the criteria for random numbers until 
the fate 1940s. At this time, computers were introduced —— 
this made the use of tables of random numbers undesirable 
mainly due to the large amount of space required for storing 
tables. In order to be able to use random digits, it was 
necessary to derive some method of generating the numbers 
within the computer itself. Two basic types of generators, 
physical and arithmetic, have been developed for internal 
computer use. Physical generators, which basically consist 
of some external device which delivers a series of random 
pulses to the computer to be converted into random numbers, 
are discussed in some detail by Scott (1967), and will not 


be elaborated on here. Arithmetic generators, which form 
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the basis for the majority of computerized random number 
generators, will be discussed in the following section of 


this chapter. 


ieee Are cnmMet lc Generators 


Arithmetic generators all have the same basic property 
= a sequence Of Gvaits is manapulaced sing a serves OT 
arithmetic operations to form a second sequence of digits. 
Part of this new sequence is retained as a new random number. 
The major advantage with this type of generator is that any 
sequence of random numbers can be reproduced exactly, if 
necessary, for checking purposes. It should be noted that 
all the numbers in a sequence of random numbers obtained 
from an arithmetic generator can be "predicted" — there- 
fore, each number in such a sequence is actually a pseudo- 


random number. 


The first arithmetic generator was developed in 1946 
by von Neumann and Metropolis. Their method, called "mid- 


square", was as follows: 


a) Define an, N-digit number x. 


b) Square this number to give a 2n digit product 
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oO) Retain the middle (nm) digats as. x and repeat 
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the process. 
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For example, consider’ =n’ =°2 ;,°"and x. = 63. Then we ob- 


tain the following: 


This generator was proven to be unsatisfactory, as the length 


of a cycle (the length of a non-repeating series) depends on 


the value of x and some initial values of x7 elead™to 
Greremely short cycles. *Forvexample ,S@consider™in =—4 55) and 
x, = 3600. The generated sequence becomes: 


een C00 PEs a= 1 OO0REX MR 2 C00mx P= 8 S00UR ee 7C00IK.h 
Another drawback is that certain starting values x lead to 
a non-repeating series of length 1 (fox example, n= 4 _, 
and os 1000) . Apart from these problems, the mid-square 


generator has been shown to not give a true uniform distribu- 


tion of digits regardlesssofvthesvalue of ox). 


Forsythe (1951) described an improved version of the 
"mid-square" generator. It is known as the "mid-product" 


generator, and is given as follows: 


a) Define two n-digit numbers x,_, , Xj - 


b) Multiply them together to obtain a 2n digit pro- 


auce. 
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Xeaqpo and repeat the process. 
For’ cxample, consider ent= 45, asa =n2344, as eooUe. “then 


we obtain the following: 


x)= 0066 x= SEE eS ah EIS: See OM Nps) 8 Aer 
This method may produce a much longer cycle than the mid- 
Square method — note that both X 5 and x, must occur 

again consecutively before the cycle can repeat. However, 


even this method has starting values x, and x, that lead 


to a non-repeating series of length 1. For example, n= 4 


0500 leads to the following series 


x O10077,, wand, ox 
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Lehmer (1949) proposed a generator that not only 
eliminated degenerating sequences, but also guaranteed a 
non-repeating cycle of maximum length. His method, called 
"multiplicative congruential", consists of forming a sequence 
of pseudo-random numbers according to the formula: 


x = ax, (Mod M) 


iad 


fOr given as and x. the maximal length for a cycle is 


ensured by the following restrictions: 
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a) x, must be relatively prime to M. 


b)y a” Ls a - primitive rooteoer b ees 9 p” Lseachac— 


COTTOr™* Me; wath? “p= ~odd=and™ a" *as® large ‘as 


possible. 
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With, 002-2 .. Moreover, forjany value.of ©M,, 


there exist values of a satisfying these conditions, 
and finally, the maximal period is the lowest common 
multiple of the periods feeaiy ee oie gene with 


respect to the prime power factors. 


Proof of the above statements is given in Hull and Dobell 
(1962). This type of generator has been widely used, and 
has passed a number of tests for randomness. Marsaglia 
(1968) found, however, that "all the pseudo-random points 
generated in the unit n-cube by a multiplicative congruen- 
tial generator will be found to lie in a relatively small 
number of parallel hyperplanes." To illustrate, suppose 
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ent points on the unit 6-cube. Marsaglia showed that if 
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all of the points were generated by a multiplicative con- 
Gruential-generator On. a, ,32-bit woxd Size computer, suchas 
the IBM 360/67, then all the points Ts would fall on less 


than 120 parallel hyperplanes |Marsaglia’s theorem states 
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that the number of hyperplanes is less than or equal to 


n ane 
(n!M) |. Tn OLchexr words, Multiplicative ‘congruential 


generators suffer from lack of multidimensional uniformity. 


Thompson (1958), Coveyou (1960), and Rotenberg (1960) 
all proposed variations to Lehmer's method, which consisted 
Of adding a constant G, to the equation) for the multiplica- 
tive congruential generator to obtain: 


x = ax. + c(Mod M) 


nee 


This variation is referred to as the mixed congruential 
method. The maximum length for a non-repeating cycle is 
ensured by the following restrictions: 
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Proot of the above statements can. be, found in Hull and DobellL 
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generators have the advantage over multiplicative congruen- 


tial generators in that a non-repeating cycle of length M 
can be attained — however, this is done at a cost of poor 


statistical behaviour (see Hull and Dobell (1964))}. 
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a) Consider the congruential generators 
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MacLaren and Marsaglia successfully developed such a generator 
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L=L*ML 
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sum (Modulo 2°?) of the randomly chosen N(J) , the 
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are the outputs of the three congruential generators, and 
the constants MK , ML, and MM are chosen to ensure 
long non-repeating cycles. Marsaglia and Bray claim to 
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Hutchinson (1966) reopened the case for using a multi- 
plicative congruential pseudo-random number generator by 
suggesting one in which the value of M (the modulus) is 
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Of generators .a) passes the besic statistical tests for 
randomness, b) has a sufficiently long non-repeating cycle 
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patterns that other congruential generators have been shown 
to Nave. Scott (1967) Gives” the venerator X44 = 24x, (Mod ey 


aS an example of Hutchinson's generator. 


Tausworthe (1965) proposed a new uniform random number 
generator, based on irreducible primitive trinomials over 
thesiszeldvor characteristic 2.. Whittlesey) (1968) “describes 
the basic properties of a Tausworthe generator, and gives the 


following algorithm for implementation purposes. 
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This algorithm will generate all possible ok non-zero 
N-bit numbers before it repeats any of these numbers for 
DrOperachowces: of. iM used withea —parnutcular SN en) For 
example, on the IBM360, where the single-precision word 
Tengen is232. b1ts o(N=3..) iM somusi be scet equals ato nse, 

7, Or 13 to obtain the maximum length of the non-repeating 
cycle. Such a generator has been tested extensively in 
recent years. Whittlesey (1968) states that a Tausworthe 
generator has no serial correlation, and speed. The same 
author (1969) states that this generator does have multi- 
dimensional uniformity.  Tootill, Robinson, and Adams (1971) 
have shown that as long as M is not too small, not too 
elose toe N72), Sor noe oreater than N72e,) this) cenerator 
also has a good runs up and down performance. However, 

the Tausworthe generator suffers from the fact that for 
many computer word lengths N , the values of M producing 
non-repeating cycles of length pore 1 are not known (only 
Zierler and Brillhart (1969) have made any attempt to tabu- 


late values of M for certain given values of N). 
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CHAPTER@2 


DESCRIPTION OF NEW SYSTEM OF GENERATORS 


Zete IRCrOduUCE ION 


In the previous chapter, the development of uniform 
pseudorandom number generators was outlined, with special 
emphasis being placed on the so-called congruential system 
of generators, and the primitive polynomial based generator 
of Tausworthe. It is apparent, however, that the optimal 
use of any congruential or Tausworthe generator is dependent 
on the word-size in bits allowed on various computers. For 


example, the multiplicative congruential generator defined 
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behaviour, while using a minimal amount of computer storage. 


2.2 Description of General System 


The basic elements for the new system are as follows: 
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2. Define new value(s) for any pointers used. 
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3. Reorder the column of A _ from which the digit was 


selected. 


There are many general procedures that could be derived 


for manipulating A in order to produce a sequence of random 


digits. Two such procedures, namely Sequential Single Arrays 


and Sequential Pair Arrays, have been analyzed in terms of 


the components given; both of these will be discussed in some 


detail in the succeeding sections. 


2.3 Sequential Single Arrays 
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while p is the value of the previously generated 
digit. The next dicit will be selected from the 
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Q will become the new value of p , to be used in 
generating the next digit. 
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initial values by the user, and the starting value of j is 
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Cases, a second column pointer Zz may be required, which 
has the same range as z ..All.of these pointers have been 


assigned values during generation of the previous digit. 
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Suppose that 


represents an adjacent column pair to be manipulated. 
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have been manipulated, the value of the selected 


Gigi el wild be: 


The new value of z will be defined in the same 
Manner as in the previous method. The algorithm 
to generate an L-digit number using the preceding 
principles is the same as that for the previous 


method, but with the following two changes: 
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The LTeordering off columns P20-! Vand 2% as now 
done using the same procedure as in the preceding 
two methods, and the value of the selected digit 

Q (as well as the new values for any pointers 
used) are defined as before. An algorithm to 
generate an L-digit number using the preceding 
principles is the same as for method a., but with 
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Each time a sequence of L digits is generated by any of 
the previously described methods, the first column of B 
can be used as an index vector to re-order the L digits. 


The macris —B ‘can then be modified as follows: 
aj interchange two random digits in the first column. 
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CHAPTER. 3 


TESTING THE NEW CLASS OF PSEUDORANDOM NUMBER GENERATORS 


3.1 Introduction - Description of Tested Generators 


A series of tests (see Appendix A) were applied to each 
of the developed methods using different values of the matrix 
A . This was done in order to enable us to compare these 
methods with other already existing methods of uniform 
pseudorandom number generation. Some preliminary testing 


was done on the following generators: 
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e) Tausworthe Generator (see Whittlesey (1968)) based 
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£) Sequential single arrays, method 1 (modified with 


Selector Matrix and revised 6 function). with on = 
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tests described in Appendix A as applied to sequences of 
L-digit numbers generated by the random number generators 
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were Generated, gandrthe average “of .the various \statistics 
calculated for each sequence were entered in the tables. For 
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Ae fe 10; 15, 20) (were calculated tor keacn sequence. 
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MC - IBM Multiplicative Congruential. 

CO - Marsaglia & Bray Composite Congruential. 

T - Tausworthe. 

M,1, M,2, M,3, M.4 - Sequential Single Arrays, methods 
1 to 3, and 1 (modified), with 3 column A matrix. 
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Table 3.3 Sample Covariances (continued) 
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The numbers in brackets represent the number of trials used 
to obtain the serial correlation when machine limitations 
limited the number of trials to less than the 100 used for 
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Due to machine limitations, the results obtained for 
a digit series length of 30000 cannot be considered as 
conclusive, as in most cases, the number of trials involved 
waS very small. In fact, whenever the number of trials was 
less than 5, the average serial correlation coefficient 
obtained was not used for the purpose of setting up Table 
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geo eCOnclusions 


the -resuttes+or =the tests as Summara zed “in Tables 3.) and 
3.2 show that both Secuential Single Array and Sequential Pair 
Array generators appear to have good statistical behaviour. 
It should be noted that a Sequential Pair Array generator 
Misinge mW “array "salves tekxes More COmMpuUte: time Lo Generate 
a random number than a Sequential Single Array generator 
using 2m arrays, without improving the statistical behaviour 
MNoLLceably. Therefore, the data in Tables 3.3°t0 3.5 compares 
already existing generators with Sequential Single Array 


generators only. 


For all series lengths of generated digits, the 
Sequential Single Array systems appear to compare very favor- 
ably with walready existing systeris wlth respect to theyserral 


correlation between consecutive pairs of overlapping digits: 
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improvements in the Sequential Single Array methods (or in 
the starting value of the matrix A used) could cause a 


further improvement in this respect. 


From observing Table 3.3, there does not appear to be a 
noticeable difference between Sequential Single Array systems, 
and already existing systems with respect to the observed 
sample covariances of generated sequences of numbers. The 
difference between the maximum positive and negative values 
for all tested lags is smaller in many cases for the Sequen- 
tial Single Array generated sequences than for sequences 
generated by already existing methods (note that zero lies 
well within the observed range of values in all cases). 

Again, improvements in the methods (or in the starting value 


of A) could cause the differences to become more noticeable. 


The same algorithms for the Sequential Single Array 
generators were run on an IBM 360/67 (32-bit word length) 
and a DEC-PDP9 (16-bit word length). As these two computers 
represent an extreme in allowable bit-word length, it would 
appear that Sequential Single Array generators should be 
able to be implemented without any changes in the algorithm 
involved, or parameter values used, on any computer that 
allows array manipulation. The same is not necessarily true 


for a Congruential or Tausworthe generator. 


The following table shows a comparison between various 


generators in terms of time and core used, using acpuale data 
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taken from testing with L = 2, and all testing subroutines 


eliminated. 


Table 3.5 9Time and Core Comparison 


Generator 


Core Used 
(bytes) 


Time refers to the time in seconds taken to generate 
100 numbers of 2 digits each. Core Used refers to the number 
Of bytes used to store the FORTRAN object code pilus all 


variables needed for generating a random number of L digits. 


It should be noted that the extremely high table values 
for the Tausworthe generator are caused by the fact that the 
FORTRAN implementation of such a generator involves simula- 


tion of the Exclusive-OR operation. The use of a language 
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such as ASSEMBLER (for which Exclusive-OR is part of the 
basic instruction set) would drastically reduce both time 

and core used. As for the other generators, both the Mult-— 
iplicative Congruential and Composite Congruential generators 
are much faster than the new generators. However, in terms 
of core used, the new systems as they now stand (particularly 
Methods 1 and 2) compare favorably with the Composite Con- 
gruential generator. It may be that more efficient program 
coding will reduce the core used for the new generators to 
such an extent that they will all save core storage when 
compared to the Composite Congruential generator, and may 
even approach the 900 bytes used for a Multiplicative Con- 


gruential generator. 


When looking at CORE USED for all of the generators, 


tWO RLaCtOrs. should be considered: 


a) For Ml, M2, and M3, as L (the number of digits used 
to form each random number) increases by 1, CORE USED will 
increase by 4 bytes (these 4 bytes being used to store the 
extra digit generated). For M4, as L increases by 1, CORE 
USED will increase by (4+4L*) bytes (the same 4 bytes as 
before, plus the extra core used to store the extra elements 


of the selector matrix B). 


b) For all the new generators, as M increases by l, 


CORE USED will be increased by the 40 bytes used to store 
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the new column of the A matrix. 


The following is an example of how CORE USED was cal- 


culated for method 3 with I=2, M=3. 


CORE USED 
BYLES EXPLANATION 
L208 FORTRAN object code 
120 A-matrix of (10x3) numbers 
4 value of L 
4 value of M 
4 column pointer 
4 special pointer 
4 column index K 
4 the final L-digit number 
8 the two generated digits, stored 
separately 


TOTAL 1880 bytes 


+tiia enw fiRet ana! wort Bo shigean cia #2 on 
eee, End Made £6 


i 


WO. TA Uy Oke 
sfc thetdo MARTHOS 
esetheugr fa ang) 10 diesen 
& Ie aulev 
Map colar 
AAPBEOS sitio 
wYaingtog ‘arcege 
i *shat antioe 
~ertetn ‘Aberin Yeats art 
hevota ,<diphb badezsse ova ads 
Ylasseqas 


as 
¥} we. 2 : 


42. 


With respect to the already reasonable behaviour of 
the new systems, further research should be done in order 
Eo determine 1) if the starting parameter walues. for a 
sequential Single Array system can be optimized for a parti- 
Gilat nunvem On. Columns Ean Mate BA, Bean aul) esi Cue ner 
increase in the number of columns in A improves the stat- 
istical behaviour of the generator enough to offset any 
increase in computer storage required. Any work done along 
these lines should take the followings tact: into, account: 
an increase in the number of columns of A causes an in- 
crease in both the amount of computer storage required, 
and the expected length of a unique series of generated 
digits. The latter part of this statement can be seen from 
the fact that for a starting value for matrix A with n 
Columns, the probabidicy that A will again take: on its 


starting value is: 


A el ete 
10 ano 
This suggests that the expected value of the unique digit 


series length is 


10 x (10!)" 
n e 
If one also takes into account the probability that 
the various pointers will also return to their starting 


values at the same time as Matrix’ A» , this expected unique 
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digit length will be increased further. 


There is another distinct advantage to the new system 
of generators. Unlike any Tausworthe or congruential gen- 
erator, the new generators can be easily adapted to the 
generation of random digits using any base arithmetic up 
to DasenlOm) the general programs for generating random digits 
using Sequential Single Array methods for any such base 
arithmetic are given in Appendix B. Only base 10 generators 
Were tested ere in order to provide divect comparisons 


with the already existing generators. 


It might be wise in the future, when analyzing the new 
GeEneracvons, tO incorporate tests of Significance in analyzing 
the results ofjspecific tests for uniformity. In this way, 
one can truly determine if a generator passes or fails certain 


eriteria. 


To summarize, in spite of an increase in computer time 
required for implementing a Sequential Single Array digitized 
generator, these new generators as tested appear to compare 


favorably with existing generators with regard to statistical 


behavrour, and, inemost Cases storage .o  Thesnewsgenerators 
are also machine-independent: this is not necessarily true 
for any congruential or Tausworthe generator. It is con- 


cluded, therefore, that Sequential Single Array generators 
do have some merit, and that any further work to be done 


shoulda concentrate On Optimizing these igenexators. 
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Appendix A 
Description of Tests for Uniform Pseudorandom 


Number Generation 


A series of tests were used in order to determine the 
randomness of a generated sequence of uniform pseudorandom 
numbers. Suppose, then, that we have a sequence of I 
numbers of length L digits each Coes ace 7 
generated by some uniform random number generator. The 
following should hold true if the generated numbers are 


fromm che tEue Um Lor ClSerLouclon-: 


a) Moments (see Gorenstein (1967) }. 
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i) The first moment, ) 


has expected value 
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iv) The sample covariance, 


ne 2L —— 
a *i*it+y / (a * 107”) ite XX ayy 


should be equal to zero for any value of | eee 


a I-y ie 
where x, = os x, (fry) maa e, 


Runs Up and Down (see Downham (1969)). 


Define r as the number of runs in the sequence of 
I numbers, and ym) as the number of runs above 
and below the mean in the sequence of I numbers. 


Then: 


a) gan Gie be = (21 - 1) 
ii) B(r'™) = $41 


Frequency Tests {see Jansson (1966)). 


i) One-dimensional. One divides the interval 
0 Sale into inves equal non-overlapping sub- 
intervals, and calculates the frequency dis- 
ip DU CLON OD Cie eel a= nays numbers with respect 


to these intervals. As the expected value of 
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the number of L-digit numbers falling into 
each interval iss 1/5) eu Sstatisticscan 
then be calculated for testing purposes by 


using 


where O; is the observed number of L-digit 
numbers falling in the ee interval, while 
ES =i] -f0r able Bathor thepurpose of test— 
ing the generators described here, 100 values 
of  y* were calculated for 100 sequences of 


I= Hite numbers, and the average of these 


Walues taken as a representative y° statis-— 


ELeewien 10" degrees of freedom, 


Two-dimensional. One divides the existence 
square Coes 0 ron) into equal non- 
overlapping Squares, and calculates the 
frequency distribution of successive pairs of 
the I numbers in these squares. These pairs 
may be overlapping [ies roa (pts i 

wee y Coat , Or non-overlapping 
Peeves seen e Bee ea eae) oeAn Oar 
of these points should show a uniform distri- 


bution over the existence square. 
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Serial Association (see Keeping (1962)). 


The’ folowing vis van salgori thm (for calculating a 


serial association coefficient for a sequence of 


iP numbers on length Ti sdigitseeach: 


ay) 


hay) 


ie.) 


iv) 


Treat the series of =I numbers as a series 


Oia Gage teed cute: 


Calculate the frequency distribution table 
Of the mumber olf times! digit value “i, jis 


followed by digit value j in the series. 


Divide all elements in this table by the 
total number of overlapping pairs of consecu- 
EiVel digits EO. ObDtainna Drobabulity sable. 
Therefore, define Mae as the probability 
that digit value 1, "cecurs , and is followed 
by digit value =, 8 and I. as the 
PrLOobavulicy ethereal eval lem 1 OCCULS = neq 


Gardless of digit value 3 


Calculate the quantity 


oie 
a a actian =| /9 
/ 


i=0 j=0 “i 5 


This quantity 1S a measure Of the association 


of consecutive digits in the series, and 
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Appendix B 


Programs 


The purpose Of this appendix is to present listings of 


program segments written in FORTRAN IV to generate an 


L-digit random number using the various models of Sequential 


Single Array generators described in Chapter 2, as well as 


listings of programs for testing sequences of numbers 


generated by such generators. »in all wcases,the.fol lowing 


variables have been defined as the result of generating some 


previous number 


(or initializing the system involved) : 


P =" the arithmetic pase value being used: 

L - the digit length of the number being generated. 
TA = ‘the matrix “A 

M = the number of columns ian A 

IP - reference pointer to an element in a column of A 
K - reference column index pointer for A 

EP28—> Sspecivalecol umn sindex poi nteneror lA 

IB - the 2xL selector matrix B 


In addition, the following two variables are assigned values 


in the course of generating a number: 


NUM — (a vector containing the &° digits of a generated 


number. 
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the final generated number, formed from the 


elements of NUM. 


For testing programs, a sequence of LM consecutively 


generated values of) KNUM)<is taken (call this sequence 


NUMB), and the following variables defined: 


EERO’ Ss a L0PG10 Smatrix, an which the et element 
will contain the number of times digit value i 
is followed by digit value j in the sequence 
Or LM i-digit numbers. 

ILN = a Vector Of Length, b> containing the (le gdigics 
of a particular element of NUMB to be used 
an derining es leRO. 

AVG, AVG2, AVG3 - the three moments. 

IRUNS - the number of runs up and down. 

IRNS - the number of runs above and below the mean. 

PHI = the serial correlation coefficient, 

COV - the sample covariance. 

VFREQ - a vector containing the frequency count used 
for caleulating <= 

CHISQ - the y* value. 


The program segments are as follows: 
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b) 


Method 1. 


No O71 


teal 


Method 2. 


N U1 W 


1 


KNUM=0 

Doe? Wale 

J=IA(K,IP) 

Th=tA (ODMR EM) uo 1) 
TOsTACh a1) 

IA(K,IP)=10 

IA (K,IL+1)=d 

NUM (N) =IQ 

IP=IQ+1 

IF(IO.NE.O) GO TO 3 

IQ=1 

DO 5 KL=1,P 

TA (K, Ri) =MOD (1A (K,Kl)+10,P) 
K=MOD (K,M) +1 

DO ldo = let 
KNUM=KNUM+NUM (JJ) *10** (L-JJ) 


KNUM=0 

DO 2 N=1,L 

J=IA(K,IP) 

NUM (N) =J 

IA (K,IP)=IA(K,P+1-IP) 

IA (K,P+1-IP)=J 

IP=J+1 

LENE 60) "GO TO=3 

J=1 

DGOvoe hie le 

TACK RL) MOD. (TALK, Kia) +d, 2) 
K=MOD (K,M) +1 

DOB ad da 
KNUM=KNUM+NUM (JJ) *10** (L-JJ) 
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Method 3. 


983 


Method 1, 


ND O1 W 


all 


96 


KNUM=0 

DO 2 N=1,L 

J=IA(K,IP) 

To=1At 1 Po ood) 

TO=2A (Kk, 1t+1) 
IF(MOD(IQ,M).NE.M-1) GO TO 3 
IP2=MOD (IP2,M) +1 

GO TO 5 

IP2=MOD (IP2+1IQ,M)+1l 
IA(K,IP)=I0 

IA (K,IL+1) =3 

IP=IQ+1 

NUM (N) =IQ 

IF(IQ.NE.0)GO TO 983 

TQ=1 

DO 6 KL=1,P 

IA (K,KL) =MOD (IA (K , KL) +10,P) 
K=MOD (K,M) +1 

Dro} alelt fafer obese, 
KNUM=KNUM+NUM (JJ) *10%** (L-JJ) 


Wlth Selector matrix. 


KNUM=0 

DOV 2a N= lt 

J=IA(K,IP) 

TL=IA (MOD (K,M)+1,J+1) 

TO=IA (K,;iL+1) 

TA(K,IP)=IQ 

TAK, Lie) =o) 

NUM(N) =IQ 

IP=IQ+1 

TE (LOCNE. 0) GO T0773 
IQ=(MOD(J,P) +1) 

DO 5 KL=1,P 

TANhy i) MOD (TAK, KL) +00, P) 
CONTINUE 

DO Wi wos 
KNUM=KNUM+NUM (IB (1,dJJ) ) *10** (L-JJ) 
1Z=I1B(2,MOD(IP,L)+1) 
TR=MOD (IP ,L)+1 

TIT=ne yay 

TE een eele,e) 
TB(L; TR) HLT 

DO 96 I=1,L 

IB) SOD (ely) NUM CS ee 
p= Le) 

DB lee) = eee a) 

ete) 1 0 
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SO) slesiSeo1 Ving» oate fora Tables. 3. l= 9.0. 


PE INUME (2) eNUMB aro 5b 


iB) IsG=-l 
SOTO 7 

LG ISG=1 

uy IRUNS=1 


DO 14 I=3,LM 
IF ( (NUMB (I) -NUMB (I-1))*ISG) 18,14,14 


18 IRUNS=IRUNS+1 
ISG=-ISG 
14 CONTINUE 
TEANUMB (1). LM/2.0:)  COeTOn32 
ISG=-1 
COMO a> 
32 ISG=+1 
sie) IRNS=1 


DO 347 =27 5M 
IF ((NUMB(I) -LM/2.0-0.5)*ISG) 35, 34, 34 
35 I RNS=IRNS+1 
PSG=—-1SG 
34 CONTINUE 
AVG=0.0 
DViG2=0 0 
AVG3=0.0 
DO 24s t= 
AVG=AVG+FLOAT (NUMB (I) ) /LM 
AVG2=AVG2+ (FLOAT (NUMB (I) ) /LM) **2 
24 AVG3=AVG3+ (FLOAT (NUMB (I) ) /LM) **3 
AVG=AVG/LM 
AVG2=AVG2/ OM 
AVG3=AVG3/UM 
DOr lela.) uM 
iL VFREQ (I) =0.0 
DO wz EM 
KT=NUMB (I) +1 
2 VFREO (KT) =VFREQ (KT) +1 
CHISQ=0.0 
DOS) fd, OM 
3 CHISOQ=CHISQ+ ( (VFREQ (I) -1) **2) 
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fe) test £0r Serval Correlation Coctivecient. 


DO: 668 a—=1710 
DO” S6eaa=—, 10 
868 TEFRO (Geo 
DO’ 50° I=i7 iM 
DO} ie R= ee 
ILN (K) = (NUMB (I) /(10**K) )*10**K 
Dil: ILN (K) =NUMB (1) -ILN (K) 
DOT SZ hors 
52 ILN (L+2-K) = (ILN (L+2-K) -ILN (L+1-K) ) /10** (L+1-K) 
DO G53 42275 
5s IFRQ (ILN (K-1)+1,ILN (K) +1) =IFRQ (ILN (K-1)+1,ILN(K)+1)+1 
LRCICEO SEM) 2 <GO TO: 50 
KR=NUMB (I+1) /10** (L-1) 
IFRQ (ILN (K-1) +1,KR+1) =IFROQ (ILN (K-1)+1,KR+1) 41 
50 CONTINUE 
FSUM=0.0 
DO m2 el leo 
DOW 12, = 10 
DAD FSUM=FSUM+IFRQ (I,J) 
POLI eI 19 0 
DOm2 Po aoe) 
PALS FMAT (TI ,J)=IFROQ (I,J) /FSUM 
DOA iG 
FPRB(I)=0.0 
DOm4 yg — 17,0 
4 FPRB (1) =FPRB(1)+FMAT (I,J) 
DOu Sno lho 
GPRB (J)=0.0 
DOm>e iy Ld 
5 GPRB (J) =GPRB (J) +FMAT (I,J) 
PHI=0.0 
DOM Gegl— 15,00 
DOnGmo— 1 0 
6 PHI = PHI+((FMAT (I,J) **2.0)/(GPRB(J) *FPRB(TI))) 
PHI=(PHI-1.0) /9.0 
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g) Test for Sample Covariance. 


DOW 550 1G=17,.20 
TE (IG Gr. LO AND. LG. dil. bo). OR (iG.Gl..U5.. 
LOR 1G lot 20)) GO TORS 50 
LP=LM-IG 
SUM1=0.0 
SUM2=0.0 
SUM3=0.0 
DO 551 KG=1,LP 
SUM1=SUM1+NUMB (KG) /FLOAT (LM) 
SUM2=SUM2+NUMB (KG+IG) /FLOAT (LM) 

551 SUM3=SUM3+NUMB (KG) *NUMB (KG+IG) /FLOAT (LM** 2) 
COV= (SUM3/LP) - ( (SUML*SUM2) / (LP**2) ) 

550 CONTINUE 
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Appendix C 


Period cmcy 


One aspect of the new random number generators that 
was not discussed previously in any detail was the actual 
periodicity of the sequences of digits generated. Due to a 
lack in availability of computer time and resources, it was 
not considered feasible to test the base 10 Sequential Single 
Array generators of Chapter 3 extensively for periodicity’. 
However, as the question of the possible effect of an in- 
crease in the number of columns of matrix A _ on the period- 
icity was raised, it was decided to perform some testing for 
periodicity on a base 4 version of one of the new Sequen- 
tial Single Array generators. Method 2 was chosen, and 


the following values of A were used: 


1A total of 60,000,000 digits were generated using the base 

10 generator referred to in Chapter 3 by M,l, without any 
repetition in the set of starting parameter values, and with- 
out any internal degenerate sequences occuring. However, due 
to the high cost involved in generating that many digits, the 


process was not carried out any further. 
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FOr each Of the starting values of Ay used, ald 
possible combinations of starting values for the column 
pointer Zz and pointer index p were tested, in order to 
determine which combination of p and Zz for a particular 
A produced a sequence of digits of maximum periodicity. 


The results were as follows: 


4-column matrices 5-column matrices 
Matrix eon Matrix abe aairg te 

Ay 24128 A, 298590 

A, 27200 A, 254730 


It would appear as though the increase in the number of 


columns of A from 4 to 5 causes a great increase in the 


figs (Bear si° to czaulea Ra lo iota toF 
Wition OnT sol Se0 55" HLis4BFH) 3 sobtentdmnos aldies Og 
3, sa620 Ne ghasesu «2 mio” bea, & see tog 
Si +76 B» Zee cs Ose c | FeOLtanianos Acs H onimxs: of 


Fy rf a . ws 
froLbGiIrsg GinmrKea bo ib 20 enteyvpes « beopboste 


_ ‘ 
- 

he! 

J q 

7 - 

- 


- ~ ra AT 
Oliot ae stow atfoser, sat 
iz 


AasSLtzanm miltioo-< al 14 Se ar foot 


vyIiothoria9 eer 240 (ange1s4 
iz nts » J ra - q mad xiladsM 7 a 
_ 7 _ a - 
OpeHeOS. > | aciac ; 
> - 


a - - - 


peyags. 
- 8 - aa 
ey de ee 
3 it ti. sasex5a2 Ine 
pen 


Spl lys 


actual periodicity — Lt 1S Expected thatsa furcher mncrease 
inthe number offcolumns- (and=the optimization of A for 
any number of columns) will cause the increase in periodicity 


to be more pronounced. 
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Appendix D 


Applications to Queueing Theory 


Cae Ineroulcrion 


One major use of Uniform pseudorandom number generation 
is in the study and simulation of queueing systems, in which 
the generation of random arrivals and departures plays an 
ImMpOreant role. 1h as Lhe purpose of this ‘section £o define 
one particular queueing model, and to describe how one can 
optimize this model by use of Uniform pseudorandom number 


generation. 


C.2 Description of Single-Channel, Single-Server, Dynamic 
Queueing Model 


A single-channel single-server dynamic queueing system 


(no feedback) can be denoted by S(T,a,6,N) pon dt Witchy 
1) T , =the processing period (or simulation period) 
is the lapse of time between the start (t=0) and 


the end (t=T) of each processing. 


2) o SOC) ete O51, ey eS VeECcCOn OL time 


dependent job arrival rates. 


Sh) Brie Ge) to =O eel LG ce VECLON .Ote tlines 


dependent job departure rates, and 
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4) N is the maximum allowable queue length. 
It should be noted that: 


I) 8(t) is andependent of a(t) , but dependent on 


the efficiency of the server. 


2) The maximum allowable queue length, N , is depen- 
dent entirely on the setup of the server. An in- 
coming job will be rejected if there are already N 
jobs waiting in the queue. N is assumed here to 


be time independent. 


When simulating such a system, the following quantities 


should be calculated: 


1) The total number of customers arriving at the system 


ARR 


2) The total number of customers rejected by the system 


3) The total number of customers entering the server 


Nie. 


4) The total number of customers departing from the 


server DEP . 


5) Thesnumberso£t timesunits: ttle fox which~thesserver 
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Thepaggregate total of the number of units of time 


all customers have waited in the queue QTM 


The following is an _algorithm.for simulating this system: 


Liss 


inzevalrze WARR, Reg, ENT? DEP, TDi, and OTM. “to be 


ZiCieOrs 


Set the index I , and the queue-length to be zero. 


Generate two normalized L-digit random numbers 
ARVL and DEPT using a sequential single array 


method. 


If DEPT is less than 8(I) , then go to statement 


LOA) Bh LSS 
If ARVL is less than a(I) , then go to statement 
iG .4.hLSe; 


If the server is busy, then go to statement 9 


Else, 


LE the’ quene—lengteh 1s mct equal to zero, then, go to 


statement .20 «, Else, 


Increase the value of IDL by one as the server is 


1dle wandegouto stacementaec. 


Increase QTM by the value of the queue-length, and 


go.to Statement. 255. 
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iz the Server is not busy, them go to statement 16 


Else, 


tncrease Lue Value (ote DEP Gby ol) 7. as a. iob 1s 


leaving the server. 


It .ARVE gis lessothan gall). , then go to statement 


15a) Bise, 


if thesqueue—Length 1s equal to zero, then go. to 


statement 8. Else, 
Go to statement 20 


Increase the values of ENT and ARR by 1 (as a 
JOD_IS arriving at the system, vandiva job. as) entering 


the server) and go to statement 9 


If ARVL is less than a(I) , then go to statement 


L523 > JbLse, 
GO to statement $8 .« 


lf the queue—-length is mot equal to Ne, then go 


to statement “21. lise, 


increase the vallie-of “REQ by 1 (as ayjob is 


rejected by the system) and go to statement 9 


Decrease the queuve-length by 1, and increase the 


VeEiemot ENT by wh (as fa joo s,entering tie 
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24. 


25's 


26. 


Dies 


server from the queue) and go to statement 9 


Increase the queue-length by 1. 


If the server is not busy, then go to statement 24 


Else, 

Increase the value of ARR by 1, and’go to state- 
ment 9 

Decrease the queue-length by 1 , increase the 
Values OL JEN Jand ARR “by 10, and go to, state— 
ment 9 . 


Increase the value of I by 1 


iP cee Varner Or 8L ssc less Cham fo, eiens go. co 


Stacement  o. Glce, 


Stop, as the simulation Is finished. 
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Ceo Optimization of the System for Given a and N 


In order to determine the optimal behaviour of the 
Sycotem Gesuribed in section C.2, 10 1s Necessary to find the 


answers to the following two questions in turn. 


a) For any system where the area under the curve de- 
fined. byeedit)s Vs 0 t “is equal) to the “area under 
the corresponding curve. defined by 6b). vs to 
what pattern of &({t) will optimize the system 


for a, given. a(t), 2 


D) S Starting with the optimal systemior “a(t)" Fand Bt) 
values determined from answering question a), how 
Will overstaffing and understaffiing change the 


behaviour of the system? 


Before attempting to answer either of these questions, one 


must first, assign arbitrary values to: 


1) LOSS — the oss an prot foreach time muni t a 


customer waits in the queue. 


id) SGAIN —, the dain gn profit for, cach time une a 


customer spends in the server. 


2 By overstaffing and understaffing, we refer to defining 
Vealwes OL 1 16-(b) euch that the curve defined by  8(t) vs 
t as either greater than or less than (in area uncer- 
neath) the curve defined by a(t) “vs t 
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ivi) IDCOST. = the upkeep cost for each time Unit the 


server 2s idle. 


LW) SVCOST’ "= the® Pabor cost for each tame unit the 


Server is busy. 


Once these variables have been given values, any curve 8 (t) 
vs t that satisfies the criteria of either questions a) 


or b) can be used to determine the behaviour of the complete 


system for given OL pie ands eN by sche tollowind: 
Lever tne tenes Set Or Values s reuCi) 6 t= nO elaine 
i ))  Perrorm a senies Of Simulations using S(T,a,6,N) , 


and calculate the average values of the parameters 


EDEN s77-OTM 
iit). Calculate the system prozic by evaluating: 


PROFIT=((GAIN-SVCOST) «( T-IDL))- (IDCOST*IDL) - (LOSS*QTM) . 


By testing various @(t) patterns in this way, one can 
eventually find a pattern that produces a maximum value of 
PROFIT for a given a it pander N sineordem CO,enswer tne 


two questions posed above. 
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